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ABSTRACT 


This  report  presents  two  methods  for  the  estimation  of  a  periodic  signal 
in  additive  noise.  Both  methods  assume  that  only  a  finite  time  sample  of  the 
signal  plus  noise  is  available  for  processing.  The  estimate  of  the  signal  is 
chosen  to  minimize  the  mean  square  error  between  the  estimate  and  the  sample 
of  signal  plus  noise.  This  is  also  a  maximum  likelihood  estimate  if  the  noise 
is  white  and  Gaussian. 

The  first  method  is  frequency  domain  analysis.  Estimates  for  the 
period  of  the  signal  and  the  complex  amplitudes  of  its  harmonics  are  derived. 

The  second  method  is  time  domain  analysis.  Estimates  for  the  period 
of  the  signal  and  for  the  waveform  of  me  period  are  derived. 

Under  the  assumption  of  white  noise  and  large  signal-to-noise  ratio, 
formulas  for  the  expected  values  and  variances  of  the  period  estimates  are 
derived.  The  estimates  for  the  period  are  found  to  be  the  same  by  both  methods. 
The  estimate  is  unbiased  and  has  a  variance  inversely  proportional  to  the 
signal-to-noise  ratio,  and  inversely  proportional  to  the  cube  of  the  number  of 
periods  in  the  given  sample.  The  expected  values  of  the  estimates  of  the  wave¬ 
form  itself  are  derived,  and  the  estimates  are  found  to  be  biased. 
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THE  ESTIMATION  OF  PERIODIC  SIGNALS  IN  NOISE 


T.G.  Kincaid*  and  H.  J.  Scudder,  III 


I.  INTRODUCTION 

In  this  report  we  consider  the  problem  of  estimaving  a  periodic  signal 
in  additive  noise.  Although  this  is  a  special  application  of  the  established 
theory  of  parameter  estimation,  (*»  2)  periodic  signals  deserve  a  detailed  treat¬ 
ment  because  of  their  frequent  occurrence.  For  example,  the  results  of  this 
study  are  useful  in  the  analysis  of  sounds  produced  by  rotating  machinery.  In 
this  application,  good  signal  estimates  can  be  used  for  source  identification  or 
malfunction  diagnosis. 

We  shall  assume  that  a  finite  time  sample  of  the  signal  plus  noise  is 
available  for  inspection.  We  shall  show  two  methods  of  forming  the  estimate 
of  the  periodic  signal  which  minimizes  the  mean  square  error  between  this 
estimate  and  the  given  time  sample  of  signal  plus  noise.  (It  is  well  known  that 
the  minimum  mean  square  error  estimator  is  also  a  maximum  likelihood  esti¬ 
mator  if  the  noise  is  white  and  Gaussian. )  The  first  method  uses  a  frequency 
domain  approach,  while  the  second  method  works  in  the  time  domain. 

In  the  frequency  domain  method,  the  first  step  is  to  estimate  the  period 
by  summing  the  power  spectrum  of  the  data  at  harmonics  of  a  sequence  of  trial 
fundamental  frequencies.  The  period  of  the  fundamental  which  gives  the  largest 
sum  is  the  estimate.  The  Fourier  coefficients  of  desired  waveform  are  then 
estimated  from  the  values  of  the  Fourier  transform  of  the  data  at  harmonics  of 
the  estimated  fundamental. 

In  the  time  domain  method,  the  period  and  the  waveform  are  estimated 
simultaneously.  The  technique  is  simply  to  choose  a  trial  period  P,  divide  the 
data  into  P  length  sections,  and  average  the  sections.  The  trial  period  which 
produces  the  average  waveform  with  the  largest  energy  is  the  estimate  of  the 
time  period,  and  the  average  waveform  is  the  estimate  of  the  true  waveform. 

Under  the  assumption  of  white  noise  and  large  signal-to-noise  ratio, 
formulas  for  appropriate  statistics  of  the  estimates  are  derived  for  both  methods. 
The  mean  and  variance  of  the  period  estimates  are  found  to  be  the  same  by  both 
methods.  The  mean  is  found  to  be  unbiased.  The  variance  is  found  to  be  inversely 
proportional  to  the  signal-to-ncise  ratio  and  inversely  proportional  to  the  cube 
of  the  number  of  periods  in  the  data  sample.  Formulas  for  the  mean  of  the 
Fourier  coefficients  are  derived  for  the  frequency  domain  estimates,  and  it  is 
found  that  these  estimates  are  biased.  In  the  time  domain  approach,  it  is 
necessary  to  assume  a  noise  distribution  to  compute  the  mean  of  the  waveform 
estimate.  On  the  assumption  of  Gaussian  noise,  this  estimate  is  shown  to  be 


*General  Electric  Company,  Heavy  Military  Electronics  Department,  Syracuse, 
New  York. 


-1- 


biased. 


In  this  report  the  following  notation  will  be  used. 

7 (t)  denotes  a  function  defined  on  the  interval  -<xt<®,  and  is  a  theoretical 

ideal, 


f(t)  denotes  a  function  which  is  equal  to  7(t)  on  a  finite  interval  T,  and 
zero  elsewhere,  and  is  the  real  signal  we  must  work  with. 

f  0  for  |  x |  >1/2 

rect(x)  =Jl/2  for  |x|  =  1/2 

1  for  |  x |  <1/2 


sinc(x)  = 


smnx 

TTX 


for  x^O 


1  for  x=  0 

(rect  and  sine  are  a  Fourier  transform  pair); 
6{t)  denotes  the  Dirac  delta  function; 


SXM  * 


0  for  t<0  and  t^x 
1  for  (Xt<x 


,  the  gate  function; 


*  between  functions  denotes  convolution; 

*  as  a  superscript  denotes  complex  conjugation; 

E  is  an  expectation  (ensemble  averaging)  operator. 

In  general,  where  small  letters  are  used  for  time  functions,  the 
corresponding  capitals  are  used  for  their  Fourier  transforms. 

Formally,  we  initiate  our  analysis  with  the  following  assumptions: 

(i)  z(t)  =  x^t)  +  n(t)  has  bandwidth  B  Hz; 

(ii)  x<)(t)  is  periodic  with  period  Po»~,  i.  e. ,  many  harmonics; 

O 

(iii)  n(t)  is  a  sample  function  of  a  zero  mean  wide  sense  stationary 

random  process  with  covariance  function  R(t-s)  and  Power 
Spectrum  P(f); 

(iv)  For  band-limited  white  noise,  R(t-s)  =  a2  sinc(2Bt)  and  P(f) 


(v)  x0(t)  w0(t-nP0)  (see  Ref.  3) 

*->n 
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where  w<)(t)  =  0;  for  t<0,  t^P0. 

(vi)  •—  J'Pc  wo(t)dt »a2 
r 0  0 

i.  e. ,  large  signal-to-noise  ratio. 

(vii)  z(t)  is  known  on  an  interval  of  length  T»P0 
i.  e. ,  many  periods. 

II.  THE  FREQUENCY  DOMAIN  APPROACH 

In  the  frequency  domain  approach  we  begin  by  characterizing  the  un¬ 
known  periodic  waveform  Xo(t)  as  a  Fourier  series.  The  optimization  procedure 
then  yields  estimators  of  the  period  P0  and  Fourier  coefficients  c<>k  of  xo(t). 

A.  The  Estimates  of  the  Period  and  Fourier  Coefficients 


Any  periodic  signal  x(t)  can  be  written  as  a  Fourier  series  as  follows 

~x(t)=[  =kSi"/p  .  (1) 

As  our  estimate  of  x0(t),  we  seek  the  particular  periodic  signal  x(t)  -  x(t)  which 
minimizes  the  mean  square  error  between  x(t)  and  the  observed  data  £(t)  on  the 
interval  -T/2  ^  t  T/2.  That  error  is  given  by 

i  rT/2 

L  =  ~  /  [z(t)  -  x(t)]2  dt  .  (2) 

J-T/2 

A 

In  Appendix  A  it  is  shown  that  the  estimates  of  the  period  P  and  the  Fourier 
coefficients  ck  of  x(t)  are  determined  as  follows. 

First  compute  Z(f),  the  Fourier  transform  of  z(t).  Then  P  is  the  value 
of  P  which  maximizes  the  expression 

i  K 

v(p)  =  f  Y,  |z(k//p)12  (3) 

k=l 

where  K  is  the  total  number  of  harmonics  that  can  be  accommodated  by  the 
bandwidth  of  z(t).  The  Fourier  coefficient  estimates  c^  are  the  values  of  c^ 
given  by 


ck  =  ak  +  ]  bk 


(4) 
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where 


ak  =  ~  Zr  (k/P)  (5) 

bk  =  |  Zi  (k/P)  .  (6) 

Zr(f)  and  Zj(f)  are,  respectively,  the  real  and  imaginary  parts  of  Z(f).  The 
estimate  x(t)  of  the  unknown  waveform  can  then  be  generated  by  the  relation 

Ht)-y  ck  ei2nk‘/P  .  (7) 

*-*k 

These  estimates  are  independent  of  the  noise  covariance  R(t-s). 

The  estimates  are  intuitively  what  we  would  expect.  The  period  is 
estimated  by  summing  the  values  of  the  power  spectrum  of  z(t)  at  the  harmonics 
of  the  fundamental  frequency  1/P.  When  plotted  as  a  function  of  P,  one  would 
expect  this  sum  to  be  greatest  when  P  is  close  to  the  true  period.  The  Fourier 
coefficients  of  x<)(t)  are  then  just  l/T  the  values  of  the  Fourier  transform  at 
harmonics  of  the  estimated  fundamental  frequency. 

B.  The  Resolution  of  the  Period  Estimate 


In  practice,  V(P)  cannot  be  tested  for  a  maximum  at  every  value  of  P. 
Usually  P0  is  known  approximately  (from  auto-correlation  or  power  spectrum 
data),  so  that  a  range  of  values  of  P  to  be  tested  can  be  established.  In  order 
to  reduce  the  computing  time,  it  is  desirable  to  test  as  few  values  of  P  as 
necessary  within  that  range.  How  far  apart  can  we  choose  the  test  values  of  P 
without  danger  of  missing  the  maximum  of  V(P)?  To  answer  this  question  we 
need  to  know  the  behavior  of  V(P)  in  the  vicinity  of  the  maximum. 


Since  the  Fourier  transform  is  a  linear  operation,  the  assumption  of 
a  large  signal-to-noise  ratio  carries  over  into  the  frequency  domain.  Thus 
the  gross  behavior  of  V(P)  in  the  vicinity  of  the  maximum  will  be  dominated  by 
the  signal,  and  we  can  study  that  behavior  adequately  by  considering  the  ideal 
case  in  which  there  is  no  noise.  In  Appendix  B  the  noise  free  behavior  of  V(P) 
is  studied  in  a  region  AP  about  P0.  For  h«P0  it  is  shown  that 


V(P0  + 


|cok | 2  sine2 


(8) 


Equation  (8)  shows  that  in  the  vicinity  of  P0  the  noise-free  V(P)  is  a  weighted 
sum  of  squared  sine  functions,  centered  on  P0.  The  widths  of  the  main  lobes 
of  the  sine  functions  are  inversely  proportional  to  k,  each  lobe  going  to  zero  at 
h  =  P§/kT.  The  "roll  off"  of  V(P)  in  the  immediate  vicinity  of  P0  is  controlled 
by  the  narrowest  of  these  squared  sine  functions.  This  is  the  term  in  Eq.  (8) 
for  which  k=K.  Thus,  in  order  not  to  miss  the  maximum,  V(P)  should  be 
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sampled  at  intervals  not  larger  than 


AP  = 

KT 


(9) 


We  shall  call  the  quantity  AP  given  by  Eq.  (9)  the  resolution  of  the  period 
estimator. 


Since  K  is  the  total  number  of  harmonics  in  ^(t),  it  is  appropriate  to 
define  trie  bandwidth  B0  of  3q>(t)  to  be  K/P0.  Then  Eq.  (9)  can  be  written 


(10) 


Equation  (10)  is  satisfying  in  its  simplicity.  It  says  that  the  resolution 
of  the  period  estimator  is  equal  to  the  true  period  divided  by  an  appropriately 
defined  time  bandwidth  product.  The  resolution  does  not  depend  upon  the  detailed 
structure  of  the  periodic  waveform. 


C.  The  Expectation  and  Variance  of  the  Period  Estimate 

Because  of  the  noise,  we  cannot  expect  the  maximum  of  V(P)  to  lie 
exactly  at  P0.  As  a  measure  of  how  close  Pis  likely  to  come  to  P0,  we  shall 
compute  the  expectation  and  the  variance  of  the  random  variable  P.  We  proceed 
as  follows.  (4) 


V(P)  is  a  well-behaved  function  of  P.  Thus  we  can  expand  V(P)  in  a 
Taylor  series  about  the  true  period  P0. 

V(P)  =  V(P0  +  h)  =  V(P0)  +  v'(P0)  h  +  V"(P0)  ~  +  ...  (ID 

We  shall  only  investigate  the  region  h  «  P0,  so  we  can  neglect  terms  of  higher 
than  second  order  in  h.  The  maximum  of  V(P)  may  be  found  by  taking  the 
derivative  of  V(P)  and  equating  it  to  zero. 

V'(P)  =  0  =  v'(P0)  +  v"(P0)h  .  (12) 

Solving  Eq.  (12)  for  h,  we  find  that  to  a  first  approximation 


£  =  Po  - 


v'CPn) 

vVo) 


U3) 


In  Appendix  C  it  is  shown  that  the  expectation  of  Pas  given  by  Eq.  (13)  is 


EP  = 


(14) 


i.  e. ,  the  estimate  is  unbiased. 
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A 

It  is  also  shown  in  Appendix  C  that  the  variance  of  P  as  given  by  Eq. 

(13) is 


Var  P  = 


(15) 


where  B2W  is  the  second  moment  of  the  spectrum  of  w0(t),  defined  by  Eq.  (C25), 
and  Ew  is  the  energy  of  wo(t),  defined  by  Eq.  (C26).  Equation  (15)  shows  that 
the  variance  of  P  decreases  as  the  inverse  cube  of  the  length  of  the  data  sample 
T,  and  that  it  decreases  as  the  inverse  of  the  signal- to -noise  ratio  Ew/o2. 

D.  The  Expectations  of  the  Fourier  Coefficient  Estimates 

The  estimates  of  the  complex  Fourier  coefficients  of  the  signal  are, 
by  Eqs.  (4),  (5),  and  (6) 


We  see  that  the  estimates  of  the  Fourier  coefficients  are  biased,  but  that  the 
bias  approaches  0  as  T-°°,  and  decreases  as  the  signal-to-noise  ratio  increases. 


III.  THE  TIME  DOMAIN  APPROACH 


We  shall  assume  in  this  analysis  that  z(t),  the  periodic  signal  plus 
noise,  is  known  on  an  interval  0<t<T.  Let  g'j'(t)  be  the  gate  function  defined  ir. 
the  Introduction.  Then 


z(t)  =  2(t)  gT(t) 


is  the  known  signal,  and 


x0(t)  =  $o(t)  g'I’(t) 
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is  the  true  periodic  waveform  obscured  by  noise 


n<)(t)  =  fy(t)  gT{t) 

on  this  interval.  In  this  approach  we  shall  find  an  estimate  w(t)  of  w0(t),  the 
waveform  of  one  period  of  3G)(t);  and  we  shall  find  an  estimate  Pof  Pq,  the  period 
of  %,(t).  These  two  estimates  give  an  estimate  x(t)  of  xo(t).  We  shall  choose 
as  our  estimate  x(t)  the  particular  x(t)  which  minimizes  the  mean  square  error 


L 


[z(t)  -  x(t)]2  dt 


(17) 


We  shall  also  determine  the  mean  and  variance  of  P,  and  the  mean  of  w(t)  under 
the  assumption  of  Gaussian  noise. 


A.  The  Estimate  of  the  Waveform  for  a  Fixed  Period 

In  this  section  we  assume  a  fixed  period  P  and  find  a  periodic  waveform 
which  minimizes  L.  To  simplify  the  analysis  we  shall  minimize  L  over  a  section 
of  data  NP  long,  rather  than  T.  The  integer  N  is  chosen  so  that  the  interval 
NP  is  as  large  as  possible  i.  e. ,  NP  4  T  <  (N+1)P.  Thus  we  shall  actually  find 
the  periodic  waveform  xp(t)  which  minimizes 


Let  Wp(t)  be  one  period  of  £p(t).  In  Appendix  D  it  is  shown  that 


Wp(t) 


z(t  +  nP)  gp(t);  NP  ^  T  <  (N  +  1)  P 


This  equation  says  that  wp(t)  is  the  average  of  N  successive  P  length  sections 
of  z(t).  This  result  is  intuitively  what  one  would  expect.  As  successive  P 
length  intervals  are  averaged  the  periodic  component  remains  constant  while 
the  noise  tends  to  average  to  zero. 


It  is  of  interest  to  examine  the  error  in  the  estimate  wp(t)  caused  by 
restricting  the  data  to  a  duration  NP.  It  is  shown  in  Appendix  D  that  the  dif¬ 
ference  between  wp(t)  given  by  Eq.  (19)  and  the  estimate  wp(t)  which  minimizes 
L  given  by  Eq.  (17)  is 


Wp(t)  -  wp(t)  = 


1 

N+l 

0 


[z(t  +  NP)  -  wp (t)j;  0  4?  t  <  (T  -  NP) 

;  (T  -  NP)  ^  t  <  T 


(20) 


This  error  is  the  difference  between  averaging  the  left  over  "tail"  of  the  data 
into  the  estimate.  Since  the  amplitude  of  error  is  inversely  proportional  to 
N+l,  our  introductory  assumption  that  T  »  P  makes  this  error  small.  By 
neglecting  this  error  we  not  only  simplify  the  analysis  but  eliminate  an  untidy 
discontinuity  in  wp(t)  at  t=T-NP. 

B.  The  Estimate  of  the  Period 

Now  that  we  have  learned  how  to  form  Wp(t)  for  an  arbitrary  P,  we 
shall  find  an  estimate  P  for  the  true  period  P0.  By  our  stated  criterion,  we 
require  P  to  be  that  P  which  minimizes  Lp  as  given  by  Eq.  (18).  It  appears 
impossible  to  get  any  satisfaction  by  setting  the  derivative  of  Lp  with  respect 
to  P  equal  to  zero.  The  P  is  tied  up  inside  the  x(t)  in  an  unknown  manner  and 
cannot  be  brought  out.  The  best  technique  appears  to  be  brute  force  computation 
of  L  as  a  function  of  P,  and  to  choose  the  P  which  minimizes  L  as  P. 

However,  it  is  not  necessary  to  actually  go  through  the  computation  of 
Lp  for  each  P  in  order  to  minimize  Lp.  It  is  sufficient  to  compute  the  energy 
of  the  estimate  Wp(t)  for  each  P.  Then  the  P  for  which  this  energy  is  a  maximum 
is  the  estimate  P  of  the  true  period  P0. 

To  show  this,  we  first  define  the  error  signal  estimate 


ep(t)  =  z(t)  -  xp(t) 

In  Appendix  E  it  is  shown  that  ep(t)  is  orthogonal  to  xp(t),  viz. , 


Xp  (t)  ep(t)  dt  =  0 


(21) 


(22) 


Now,  our  objective  was  to  choose  P  to  minimize  the  mean  square  error  signal, 
Lp. 


Since  the  first  term  is  very  nearly  a  constant  for  a  small  rang^e  of  P,  Lp  will 
be  a  minimum  when  the  second  term  is  a  maximum.  (5)  Thus  P  is  the  value  of 
P  which  maximizes  the  expression 
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V(P) 


rNP  r] 

-l  *1 


The  choice  of  notation  is  deliberate.  In  Appendix  F  it  is  shown  that  Eqs.  (24) 
and  (3)  are  equivalent  definitions  of  V(P). 

Since  the  estimator  of  the  period  is  the  same  for  both  the  time  and 
frequency  domain  approaches,  the  resolution  of  the  period  estimator  is  given 
by  Eq.  (9)  cr  Eq.  (10).  Similarly,  the  expectation  and  variance  of  the  estimate 
are  given  by  Eqs.  (14)  and  (15),  respectively. 

C.  The  Expected  Value  of  the  Waveform  Estimate 

In  this  section  we  shall  find  the  expected  value  of  w(t),  the  estimate  of 
the  waveform  of  one  period.  Unfortunately,  it  does  not  seem  possible  to  do  this 
without  assuming  a  probability  distribution  for  the  noise.  Under  the  assumption 
of  a  Gaussirn  distribution,  it  is  shown  in  Appendix  G  that 


E  w(t)  = 


x(t  +  n  P0)  s(t);  t  ^  0 


;  t  <  0 


where 


s(t)  =  1 


9 

and  erf(x)  is  the  error  function,  ai,  is  the  variance  of  P.  The  function  s(t)  is 
plotted  in  Fig.  1.  The  function  s(tf  is  a  sort  of  "soft  gate,  "  which  approaches 
the  gate  function  gp  (t)  as  the  noise  variance  a2  approaches  zero.  This  "soft¬ 
ness"  of  the  gate  reflects  the  uncertainty  in  the  period  of  x(t). 

$({) 

I  =  too 


_ yp0s  io 

— *pV  1 


Fig.  1  The  expectation 
of  the  gate  function  for 
various  ratios  of  period 
variance  to  true  period. 
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IV.  CONCLUSION 


Two  methods  of  finding  the  least  mean  square  estimate  of  a  periodic 
function  in  noise  have  been  derived.  These  have  been  labeled  the  frequency 
domain  and  time  domain  methods.  Under  appropriate  assumptions  about  the 
noise,  formulas  for  the  expectation  and  variance  of  the  period  estimate,  and  the 
expectation  of  the  waveform  estimates  have  been  found.  The  variances  of  the 
waveform  estimates  were  not  found. 

Of  the  two  methods,  the  authors  have  found  the  time  domain  method 
more  satisfactory  if  it  is  not  necessary  to  compute  the  Fourier  transform  for 
other  reasons.  If  n  data  waveform  samples  are  available,  and  it  is  desired  to 
search  m  periods  for  a  maximum,  then  the  time  required  is  proportional  to  mn. 
The  time  domain  method  has  the  additional  feature  that  the  estimates  can  be 
made  in  one  pass  through  the  data.  This  is  important  if  the  data  are  stored  on 
tape,  or  some  other  slow  memory. 

When  the  Fourier  transform  of  the  data  has  to  be  computed  anyway  for 
reasons  other  than  waveform  estimation,  the  frequency  domain  method  is 
advantageous.  Taking  the  Fourier  transform  of  n  data  points  requires  a  time 
proportional  to  n  log  n,  assuming  the  fast  Fourier  transform  is  used.  Then, 
once  the  power  spectrum  is  formed,  it  can  be  searched  exhaustively  for  evidence 
of  periodicities.  This  search  requires  mh  operations,  where  h  is  the  average 
number  of  harmonics  summed  on  each  period  tested.  This  is  essentially  a 
two- pass  system.  On  the  first  pass  the  data  are  Fourier  transformed,  and  cn 
the  second  pass  the  power  spectrum  is  searched. 

It  is  worth  noting  that  the  period  estimator  also  serves  as  a  suboptimal 
detector  of  periodic  waveforms.  When  there  are  a  large  number  of  strong 
harmonics,  it  is  somewhat  more  sensitive  than  a  spectrogram,  which  does  not 
sum  the  harmonics.  The  optimum  letector  takes  advantage  of  the  waveform 
estimate  in  an  estimator- correlator  configuration.  This  will  be  shown  by  the 
authors  in  a  forthcoming  report. 


APPENDIX  A 


We  wish  to  verify  that  P  is  the  value  of  P  which  maximizes  V(P)  given 
by  Eq.  (3),  and  that  ck  is  the  value  of  ck  given  by  Eq.  {4)  of  the  main  text. 


From  Eqs.  (1)  and  (2)  we  have  that 


L  .11 

-  T/2 

T  J. 

-T/2 

Thus, 

T  — 

1  f 

Li  - 

T  l 

|V(t)-2z(t)^r  ck  e^2nkl/P  +£  £  ckc-te 


j2TT(k+4,)t/P 


z2(t)dt  Y  ck  sine  l(k+t)  T/Pj 


-iL'.f, 


r(t)ei2nkt/P  dt 


dt  . 
(Al) 


If  we  define 


Z(k/P)  =  f  z(t)  e“32TTkt/p  dt, 

'-T/2 

the  Fourier  transform  of  the  input  z(t)  at  frequency  k/P,  and  we  use  the 
assumption  T»P,  then 

,  r  T/2  9  _  * 

L  s  -  /  z2(t)dt  +  \  |ck|2--7  cR  Z^k/P) 

1  J-T  lO  *—•  lr  1  ‘—\r 


since  only  the  diagonal  terms  of  the  double  sum  make  major  contributions. 
Lee  ck  =  ak  +  jbk,  and  note  that  bk  =  b  _k,  ak  =  a_k. 


Then 
L  5 


T  /  2 

i  f  z2(t)dt  -  Y  (a2  +  b2  -  |  y  [akZ  (k/P)-bkZi(k/P)] 
“'-T/2  £~k 


Minimizing  with  respect  to  the  ak  and  bk  yields 


ak  *  — <k/p> 


bR  »  —(k/P)  . 


Substituting,  we  find 

/T/2 

z2(t)dt  +73lY  |  Z(k/P)  j2  -  Ly  |Z(k/P) 

T  /  9  1  *->\r  1 
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or 


,  c  T/2  ,  _ 

i  z2(t)dt  -  ^  Y  |Z(k/P)|2  .  (A6) 

T  J-T/2  1 

Since  the  first  term  on  the  right  of  Eq.  (A6)  is  a  constant,  L  is 
minimized  when  the  second  term  on  the  right  is  maximized.  Since  z(t)  has  some 
finite  practical  bandwidth,  it  is  not  necessary  to  sum  over  all  k,  but  rather  to 
some  upper  limit  K.  Thus  the  value  of  P  which  maximizes  V(P)  is  given  by 
Eq.  (3).  Equations  (A 5)  give  the  values  of  and  b^  which  minimize  L  for  any 
choice  of  P.  Thus  L  is  minimized  with  respect  to  both  P  and  c^  when  c^  =  c^, 
as  given  by  Eq.  (4). 
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APPENDIX  B 


We  wish  to  verify  that  the  behavior  of  V(P)  in  the  vicinity  of  the  maxi¬ 
mum  '  i  given  by  Eq.  (8)  in  the  main  text. 

Since  we  are  assuming  the  ideal  case  of  no  noise,  z(t>  =  Xo(t).  The 
signal  is  available  between  -T/2  and  T/2,  so 

XoU)  =  rect  j^;|  £  c0n  e^2nnt^P°  .  (Bl) 

The  Fourier  transform  of  z(t)  is 


Z(f)-T^  con  sincTjf  - 


(B2) 


Substituting  Eq.  (B2)  into  Eq.  (3)  for  V(P)  gives 
K 


vtp>.  e  y  y  com c™  sincT(p''^j  s:ncT|,''~ 


t 


n 

Po 


(B3) 


If  we  set  P  =  P0  +  h  and  just  consider  the  region  where  h  <<  P0  then  Eq. 
(B3)  can  be  written 


K 


V(P0  +  h)  = 


IX  I.  com  c*nSinc  ^  k'm'l^ 

c=l  m  n  v 


T 

sine  — 
0 


.  kh 
k-n-- 


(B4) 


For  T/P0  >>  1  only  the  terms  for  which  m  =  k  and  n  =  k  will  contribute 
significantly  to  the  sum  on  k.  Setting  the  terms  for  which  m,  n/k  equal  to 
zero,  Eq.  (B4)  becomes  Eq.  (8). 
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APPENDIX  C 


We  wish  to  verify  that  the  expectation  of  P  is  given  by  Eq.  (14),  and 
that  the  variance  is  given  by  Eq.  (15).  Since  the  expression  for  V(P)  involves 
the  Fourier  transform  Z(f)  of  z(t),  it  is  a  necessary  preliminary  to  investigate 
the  Z(f). 


From  assumption  (i)  of  the  Introduction, 

Z(f)  =  X0  (f)  +  N(f)  (Cl) 

where  X0(f)  and  N(f)  are  the  Fourier  transforms  of  xo(t)  and  n(t),  respectively. 
For  the  transform  of  the  signal  alone  we  have 


X0(f)  =  J  ^  x<)(t)  rect(t/T)  e  dt 


J’-°L  [wo<‘>*^  Mt-nP)]  rect(t/T)  e~j2Tmdt 

[Wo(0  •  ££  a(f-f-|]  «Tsi»cTf 
n 

fZ  Nil  8‘"cT(f  ■ nl 

n  ' 


(C2) 


We  now  consider  the  behavior  of  the  noise  process  N(f),  which  is  a 
complex  random  variable  for  each  f.  Its  mean  is 

EN(f)  =  E  f  n(t)  e"j2TTU  dt  =  f  ^  Eff(t)  e"i2nft  =  0  . 

J-T/2 


(C3) 


The  average  noise  energy  density  is  given  by  E|N(f)|2,  and  is  determined  as 
follows: 


E  |  N(f)  | 2  =  E 


J^Eftt)  rect(t/T)  e  dt 


12 


-  E  f0  f  n(t)n(s)  rect(t/T)  rect(s/T)  e  ^2n^t  dtds 

J— ac  oo 

T/2  r  T/2 
T/2  -'-t/: 


-L 


f  i  R(t-s)  e"32TTf(t"s)  dtds 
•'-T/2 


(C4) 
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For  band- limited  white  noise 


E  j N(f) | 2  =  J*^  rect(t/T)  rect(s/T)  a*  sinc[2B(t-s)]  e  ^TT^‘  dtds 

=  a2  J^ds  rect(s/T)  e32n*w  [sinc(2Bs)*rect<s/T)  e+^2rrfS]  .  (C5) 


Taking  the  Fourier  transform  of  the  convolution,  and  integrating  over  s  first, 
we  get 


E|N(f)|2  =  2_3L  J®^  8inc2T(f-g)  rectj^j  dg  .  (C6) 

A  plot  of  l/TE|N(f)|  2,  the  average  power  density  of  the  time  limited  noise,  is 
shown  in  Fig.  C-l  for  several  different  BT  products.  Note  that  as  BT-*  ®  this 
function  approaches  the  "rect,  "  and  for  BT  >  10  a  "rect"  is  a  good  approxima¬ 
tion  to  the  actual  power  density.  This  can  be  seen  directly  from  Eq.  (C6)  by 
noting  that  compared  to  the  "rect"  function  the  squared  "sine"  function  is 
becoming  very  narrow,  i.  e. ,  approaching  an  impulse  as  BT-*®.  Since  this 
"impulse"  has  area  1/T,  we  can  use  the  approximation 

E | N(f) | 2  =  ^  rectjX)  .  (C7> 


With  these  preliminaries  taken  care  of,  we  can  now  turn  to  the  verification  of 
Eqs.  (14)  and  (15). 


f 


Fig.  C-l  Average  power  density  of  band-limited 
white  noise  for  various  BT  products. 
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From  Eqs.  (12)  and  (13)  we  note  that  v’(Pq)  and  v"(Po)  are  needed. 


v'(p0) 


A 

L,.  dP 
k 


For  large  signal-to-noise  ratios,  the  last  term  is  much  larger  than  the  middle, 
and  we  will  drop  the  middle  term.  Since  Pc  is  the  true  period,  |  X(f)|2  will 
have  a  maximum  at  each  frequency  k/P0,  so  the  first  term  is  zero.  We  are 
left  with 


V,(Po)  ~  Y*  dP  ^ Re[X(k/P0)N(k/P0)l  .  (C9) 

Also,  for  large  signal-to-noise  ratios,  we  may  neglect  the  contribution  due 
to  the  noise,  compared  with  that  of  the  signal  in  calculating 

•  <cio> 
1C 


We  work  on  v'(P)  first.  Let 

A(f)  =  2Re(X(f)N(f)]  =  X*(f)N(f)  +  N*(f)X(f)  (Cll) 

so 

~A(f)  =  A*(f)  =  X*(f)N*(f)  +  X*'(f)N(f) 

+  X(f)N*'(f)  +  X*(f)N*(f)  .  (Cl  2) 

Now  both  the  real  and  imaginary  parts  of  X(f)  are  at  a  maximum  at  the 
frequencies  k/P0,  so  X'(k/P0)  =  0  and 

A'(k/P0)  =  X*(k/Po)N,(k/P0)  +  X(k/P0)N*'(k/P0) 

v'<P0>*£  -~-A(k/Po)=£  jpr  A'(k/P0>  .  (C13) 

A  (f)  is  a  random  variable,  linearly  related  to  N* (f) 

N(f)  =  rect(t/T)n(t)  e  dt 


! 
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so 

N  (f)  =  rect(t/T)n(t)  (-j2nt)  e  ^'2TT^  dt 

Thus 

EN'(f)  =  J_m  rect(t/T)En(t)(-j2TTt)  e  ^2TT**  dt  =  0 


Hence  EV(P0)  =  0  .  (C14) 

Since  v"(P0)  as  given  by  Eq.  (CIO)  is  not  a  random  variable,  Eq.  (C14)  shows 
that  the  expectation  of  the  second  term  on  the  right  of  Eq.  (13)  is  zero.  This 
verifies  Eq.  (14). 

A 

Turning  now  to  the  variance  of  P,  we  see  from  Eq.  (13)  that  it  depends 

upon 


Var  V'(P0)  =  E[v'(P0)]2  =  £  ^ 


k  l 


-k 

pT 


-t 

p r 
r  0 


a'(tJaU) 


(C15) 


Now 

E[A'(f)A'(g)]  =  Xq  (f)Xo  (g)E[N,(f)N,(g)]  +  X0(f)  X0(g)E[N»*(f)N'*(g)] 

+  Xq  (f)Xo(g)E(N'(f)N*(g)]  +  X0(f)Xo  (g)E[N,*(f)N,(g)]  .(C16) 

We  will  evaluate  one  of  these  terms 


E[N'(f)N'(g)]=  E  rect(t/T)(-j2rrt)n(t)  e"j2TTft  dt  • 

J_ro  rect(s/T)  (-j2rrs)ff(s)  e  •'2n®s  ds 

=  J  .a,  Loo  (“j2TTt)(“j2rrs)  rect(t/T)  rect(s/T)R(t-s)  e  i2TT(ft+gs) 

=  J  (~j2nt)rect(t/T)  e  ^2rr^  dt  •  J’_ro  (-j2ns)rect(s/T)  e  ^2rr^S  R(t-s)ds 

=  j“co(-j2nt)rect(t/T)  e'‘i2TTft  dt  •  T2sinc'T(h+g)P(h)  e^2nht  dh 

=  J^T2sinc'T(h+g)P(h)dh  •  (-j2nt)rect(t/T)  e"32TT<h’f)t  dt 

=  T2sinc'T(h+g)  •  T2sinc*T(h-f)P(h)dh 


dtds 
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Substituting  for  P(h)  from  assumption  (iv)  in  the  Introduction, 

E[N’(f)N'(g)]  -  ^  T3  sinc"T(gff)  rect(g/2B)  rect(f/2B)  (C17) 

because  the  sine  function  behaves  similar  to  a  6  function  for  large  BT.  Since 
the  above  is  real,  E[N*(f)N*(g)]  =  EJN'ffjN^g)].  In  a  similar  manner. 


E[N'(f)N*'(g)]  =  EJN^tfm'tg)]  ^  ~ 
Thus 


T3sinc*T(f-g)  rect(g/2B)  rect(f/2B)  . 

(C18) 


E[A'(f)A'(g)]  a:  {[X^fjX^ (g)  +  X0(f)X0(g)]  sinc"T(f+g)  . 

+  [X0(f)x£<g)  +  X* (f)X0(g)]  sincnT(f-g)} 


rect(f/2B)  rect(g/2B)  .  (C19) 

For  f=g, 

2  0 

E[A*(f)]2  =  2 1  X0(f >  | 2  sinc"(0)  rect(f/2B) 

+  2Re[Xg(f)2fll  sinc"(2Tf)  rect(f/2B)  .  (C20) 


For  f  =  -g  we  get  the  same  expression,  so 


k|  t  \l\ 


where  we  have  summed  over  both  diagonals  separately.  We  may  neglect  the 
last  two  sets  of  terms,  in  comparison  with  the  first  since  they  are  down  by  a 
factor  of  at  least  P0/T.  Thus 


Var  V'(P) 


2-rr2  Tf.  a2 
3  B  Pj 


l  * 

|k|<BP0 


(C21 ) 
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The  denominator  of  the  second  term  on  the  right  of  Eq.  (13)  is 


vw  i2  iPo 


which  can  be  shown  to  be 


■(p„)  -r  £  ixoifj  i2”  £ 


2T2tt2 


n  *  ~ 


(C22) 


Thus  the  variance  of  the  estimator  of  the  period  P  is 


Var  *  *  VarV'(Po) 


2TT2  T3  a2 


3  B  P| 


k2lx„|JM  l! 


°  I  k  !  <BPn 


4TV 

9P§ 


k2!X0^-|2 


£  <BPn 


_  3P£  _J _ 

-  ^  2bt  y  k2|Xo(|_, 

Jk  f<BP0  P° 

Now  |X0(~)  |2  =  (T2/P§)  |  W0(k/P0)|2  and  since 

4“  Z  (k2/Po)  |W0(k/P0)  |2  *  JB  ^IWoCf)  |2  df 

P°  I  k  I  <BP  'B 


(C23) 


(C24) 


Note  that  the  integral  is  the  second  moment  of  the  power  spectrum  of 
w(t)  about  f=0.  We  could  appropriately  define  the  bandwidth  Bw  of  w(t)  by  the 
positive  square  root  of 


B2  =  - 
w  E 


r-  J*  f2  I  W(f)  |2  df 


where 


Ew  =  J.p  |W(f)  |2  df  =  JP  w2(t)  dt 


(C25) 


(C26) 


Then  the  integral  on  the  right  of  Eq.  (C24)  is  E,y.  Substituting  this  result 
back  into  Eq.  (C23)  gives  Eq.  (15). 
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APPENDIX  D 


We  shall  verify  Eqs.  (19)  and  (20).  We  shall  first  verify  Eq.  (19)  by 
finding  the  periodic  waveform  xp(t)  which  minimizes  Lp  given  by  Eq.  (18).  Let 

(N-l)P 

XP(t)  =  Y  wp(t-nP) 
m=0 

where  wp(t)  =  0;  t<0,  t  >P 
Now  let  Tj(t)  be  a  function  like  wp(t) 
i.  e. ,  q(t)  =  0;  t  <  0,  t  $.  P 


Then  let 


wp(t)  =  wp(t)  +  en(t) 

where  c  is  a  real  variable.  Since  wp(t)  minimizes  Lpf  it  is  necessary  that 

r  N 

LP(e)  =  Np“  If  z(t)  ~  £  Wp(t-nP)  -  eri(t-nP) 

L  n=0 


1  2 


dt 


be  stationary  with  respect  to  an  incremental  change  in  e  when  e  =  0.  Thus 


.  Lp(e) 

°  -  -2  IT" 

_ 1__  pNP 

=  NP  JO 


c  =  0 


N 


z(t)  -  £  wp (t-mP) 

m=0 


N 

Y  Ti(t-np)  dt 
n=0 


Now  if  we  multiply  out  the  integrand  and  use  the  fact  that  wp(t-mP)  and  rj(t-nP) 
are  in  disjoint  P  length  intervals  when  m^n,  the  right-hand  side  becomes 


N  N 

NP"  jf  '  Z  z(t)ri(t-nP)  -  Y  wp(t-nP)  ri(t-nP) 


n=0 


n=0 


dt 


N 


i  jf  lz<*>  “  Wp(t-nP]  Ti(t-nP)  dt 
n=0 


Now  make  the  substitution  of  variables 


s  =  t-nP 
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Then  the  integral  becomes 


N 

KP  £  InP  nP  [z(s+nP)  -  wp(s)]  n(s)ds 
n=0 


Since  r|(s)  is  zero  except  in  the  interval  0^  s  <  P,  this  expression  can  be 
v/ritten 


N-l 


^pT  Z  Jq  [z(s4nP)-wp(s)]  Ti(s)ds 


n=0 


Taking  the  summation  back  inside  the  first  integral,  and  making  use 
of  the  gate  function  notation,  the  expression  becomes 


i-  r1 

\tp  Jn 


rp 

NP  ^0 


rN-1 

Z  z(tf-nP)  gp(t)  -  Nwp(t) 
n=0 


h(t)  dt 


clear  that  if  this  expression  is  to  be  zero  for  all  possible  r;(t),  then  it  is  required 
that 

i  N;1 

%><*>  =  [}  L  z(t+nP >  gp(t) 

n=0 

Or,  stated  in  another  form 

x  N-‘ 


wp(t)  =  jj  Y  z(t+nP) 
n=0 


; 0^  t  <  P 

which  is  Eq.  (19).  This  is  indeed  a  minimum  since 


(Dl) 


S2Lp(e) 

71 - 


>0  for  ri(t)  /  0 


e>0 


In  order  to  verify  Eq.  (20)  it  is  necessary  to  first  repeat  the  above 
minimization  process  over  the  data  intervals  T,  rather  than  NP.  The  result  is 


wp(t) 


x 

N+l 


N 


N 


V 
/ , 

z(t+nP) 

n=0 

N-l 

r. 

z(t+nP) 

n=0 
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;  0  s£  t  <  A 


;  A  ^  t  <  P 


(D2) 


This  equation  shows  that  Wp(t)  differs  from  Wp(t)  only  in  that  the  "left-over" 
portion  of  z(t)  between  NP  and  T  has  been  averaged  in  also.  The  error  between 
the  two  estimates  is  given  by 


Wp(t)  “  Wp(t) 

N  N-l 

=  im  £  z(t+nP)  '  n  £  z(t+nP) 

n=0  n=0 

■  z(;+np)  +  nTI  -  5)  £„  z(t+nP> 

I  n-0 
N 

-  jfH  [^‘+NP>  -  n  ?  z(t'“P)] 

■  Wt+NP)  -  wp(t>] 


;  0  ^  t  <  (T-NP) 


;  0  4  t  <  (T-NP) 


;  0  .<  t  <  (T-NP) 


;  0<t<  (T-NP) 


and  zero  elsewhere.  This  is  Eq.  (20). 
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APPENDIX  E 


Equation  (22)  is  to  be  verified.  First  we  note  that 


J^  xp(t)  ep(t)dt  =  xp(t)  [z(t)  -  xp(t»  dt 


pNP  « 


xp(t)z(t)dt  -  jJ*P  Xp(t)  dt 


Second,  Xp(t)  can  be  expanded  as  follows: 


N-l 

Xp(t)  =  £  wp(t-mP) 


N-l  N-l 

N  Z!  I!  z{t+(n-m)P]  gp(t-mP) 


m=0  n=0 


Now  the  two  terms  on  the  right  of  Eq.  (EJ)  will  be  shown  to  cancel.  Expanding 
the  first  of  these  with  the  aid  of  Eq.  (E2)  gives 

pNP  1  N_1  N'1  NP 

Jq  z(t)xp(t)  dt  =  -  £  Z  Jo  z<t)zft+(n“m)p]  gp(t-mP)  dt 

m=0  n=0 

N-l  N-J 

=  N  Z  Z  C  z(t)z[t+(n-m)P]  dt  .  (E3) 

m=0  n=0 

The  second  term  on  the  right  of  Eq.  (E2)  gives 


Jo  xP(t)dt 


=  N  JQ  w^(t)  dt 


=  N  J0  N  Z  z(t+nP)  gp(t)  dt 

L  n=® 

1  N~X  N_1  pP 

=  N  Z  Z  Jo  z<t+mP)  z(t+nP)  dt 

m=0  n=0 


Substituting  s  =  t+mP 


N-l  N-l 


=  I  Z  Z  JmP+P  z<s>z[s+(n~m>p]  ds  •  <E4> 

m=0  n=0 

Since  the  right  side  of  Eq.  (E3)  is  identical  to  the  right  side  of  Eq.  (E4),  therefore 
the  right  side  of  Eq.  (El)  is  zero.  Thus,  Eq.  (22)  is  verified. 
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APPENDIX  F 


We  wish  to  show  that  Eqs.  (24)  and  (3)  are  equivalent  definitions  of 
V(P).  From  Eq.  (24), 


V(P)  =  j”P  xf>(*)dt  =  J_*  |Xp(f>ta  df 


The  second  equality  is  from  the  Parseval  theorem.  Now  by  definition 

N-l 

xp(t)  =  Wp(t-nP) 

n=0 

=  gNP(t)  [wpft)*)^  6(t-nP)] 

n 

=  gNp(t)(  ^[*(0*7  6(t+mP)l  6(t-nP)j  . 

\  l  ^m  J  'n  I 

Taking  the  Fourier  transform  of  both  sides  gives 

Xp(f)  =  NP  sinc(NPf)  e"jTTNPf*  j[z(f)]T  |  6(f  ’  ||] 

*p  sin^(Pf)  e^j^^lf-^l 

=  sinc(NPf)  e"jnNPf*  zj|)  sine  p[f  -  |j  e'jTTlf  "  p)  ]T  6  (f  -  ^) 

=  sinc(Tf)e'jTTNPf*^  ^  Zj|)  sincPj—-)  e“jnP("p“l 

k 

a  ]T  Zj|)  sine  NP  (f "  pj  e“jTTNP(f  “  p)  .  (F2! 


-jnNPf* 


Then 


4«>  -  T  Y  Z('!)z*(|j  sincNP(f-|)  sincNP(f-| 


V(P)  ■  S I  XP  (f)  I2  «  ■  nF  Yi  |zlll 


which  was  to  be  shown. 
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APPENDIX  G 

t  A 

In  tnis  Appendix  we  shall  verify  Eq.  (25).  Once  a  Phas  been  chosen 
as  the  period  estimate,  the  waveform  estimate  is 

.  N-l 

w(t)  »  -  £  z(t+nP)  gp(t)  .  (Gl) 

n*0 


Letting  P  «  P0+fi,  we  can  expand  z  in  a  Taylor  series  about  Pq. 


1  ^ 

w(t)  ■  -  Y.  (®(t+nP0)  +  z’(t+nP0)nh  +  . . .  ]  g£(t) 


For  large  signal-to-noise  ratio,  we  can  assume  that  nh«  1,  in  which  case 


N-l 

w(t)  “  -  Y  [*(t+nP0)]  gf»(t) 

n=0 


The  expectation  of  w(t)  is  then 


N-l 

Ew(t)  =  ±  Y  x<t+nPo)  Eg£(t)  +  E[n(t+nP0)  gp(t)]  .  (G2) 

n=0 


Now  we  shall  evaluate  each  of  the  expectations  on  the  right  separately. 


Egp(t)  =  £,  ep(t)  p(p)dP , 

A 

where  p(P)  is  the  density  function  of  the  period  estimate.  Unfortunately,  it  is 
necessary  to  have  a  probability  distribution  for  P  to  evaluate  Eg£(t).  Therefore, 
we  shall  proceed  with  the  assumption  that  the  noise  is  Gaussian.  Let 

A  A 

P  =  P0  +  h 

From  Eq,  (13)  we  have  that 


V'<P°> 

V"(P0) 


By  our  approximations  Vm(Pq)  is  given  by  Eq.  (C9)  and  is  not  a  function  of  the 
noise.  However,  Eqs.  (C8)  and  (C7)  show  V'(P'o)  to  be  a  linear  function  of  the 
noise.  Thus  h  is  (to  a  first  approximation)  a  Gaussian  random  variable,  with 
zero  mean  and  variance  a|>  given  by  Eq.  (15).  Then 

p(P)  =  cPp0,ap  (P> 
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1  -  §  (i^a.1  for  t£0 

1  °p  I 

0  for  t<0 

where  t  is  the  Gaussian  distribution  function.  This  allows  us  to  evaluate  the 
first  term  of  Eq.  (G2).  The  second  term  of  Eq.  (G2)  is 

E[n(tfnP0)  gp(t)] 

A 

which  is  essentially  the  average  value  of  the  noise  over  an  interval  P.  Since  the 
expected  value  of  the  noise  is  zero,  it  is  reasonable  to  expect  this  term  to  be 
nearly  zero.  Hence 

for  t>C 


for  t<0 
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5.  It  can  be  shown  that  the  change  in  the  first  term  with  P  becomes  important 
when  z2(NP)  exceeds  the  mean  square  value  of  z(t)  by  a  factor  of  N2  or 
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This  report  presents  two  methods  for  the  estimation  of  a  periodic  signal  m 
additive  noise.  Both  methods  assume  that  only  a  finite  time  sample  of  the  signal  plus 
noise  is  available  for  processing.  The  estimate  of  the  signal  is  chosen  to  minimize 
the  mean  square  error  between  the  estimate  and  the  sample  of  signal  plus  noise. 

This  is  also  a  maximum  likelihood  estimate  if  the  noise  is  white  and  Gaussian. 

The  first  method  is  frequency  domain  analysis.  Estimates  for  the  period  of  the 
signal  and  the  complex  amplitudes  of  its  harmonics  are  derived. 

The  second  method  is  time  domain  analysis.  Estimates  for  the  period  of  the 
signal  and  for  the  waveform  of  one  period  are  derived. 

Under  the  assumption  of  white  noise  and  large  signal-to-noise  ratio,  formulas 
for  the  expected  values  and  variances  of  the  period  estimates  are  derived.  The 
estimates  for  the  period  are  found  to  be  the  same  by  both  methods.  The  estimate  is 
unbiased  and  has  a  variance  inversely  proportional  to  the  signal-to-noise  ratio,  and 
inversely  proportional  to  the  cube  of  the  number  of  periods  in  the  given  sample.  The 
expected  values  at  the  estimates  of  the  waveform  itself  are  derived,  and  the 
estimates  are  found  to  be  biased. 
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